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Ph I Abstract 

. . In this article we consider the backwards diffusion problem for a 

"^ I causal diffusion model developed in [25]. Here causality means that 

2 ■ the speed of propagation of the concentration is finite. For the investi- 

gation of this inverse problem, we derive an analytic representation of 
the Green function of the causal diffusion model in the k — t— domain 
CN ■ (wave vector-time domain). We perform a theoretical and numerical 

comparison between standard (and noncausal) diffusion with our diffu- 
V/~j . sion model in the k — i— domain and in the x — t— domain. Moreover, 

C7^ i we prove that the backwards diffusion problem of the causal diffu- 

sion model is ill-posed, but not exponentially ill-posed. In contrast 
to the classical backwards diffusion problem, the forward operator of 
the causal direct problem is not compact if the space dimension is 1. 
The paper is concluded with numerical simulations of the backwards 
diffusion problem via the Landweber method. 

X: 

.<^! 1 Introduction 

The standard model of the backwards diffusion problem is a case example 
of an exponentially ill-posed problem and it is related to several problems in 
tomography and image processing. For example, such problems have been 
studied in the articles [9l|8l[71[23l[29l[l3l|2l[l8l[33l[T]and books [211 [TOl [SI 
E^lEaiSSlEHllIlEniEIlEIlto name but a few. 

This article starts over the backwards diffusion problem by replacing the 
standard diffusion equation by a causal diffusion model. Recently, such a 
model has been developed and studied in [25]. By causality we understand 
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that a characteristic feature of a process hke an interface or a front must 
propagate with a finite speed c. This means that if a "point concentration" 
is added to a solution in point x = 0, then the concentration is zero outside 
the ball Bct{0) after the time period T. Here it is not relevant whether 
this interface or front is visible. Although a causal behavior is naturally de- 
manded for problems involving hyperbolic equations, it is usually disregarded 
for problems involving parabolic equations. It seem to the author that the 
modeling of causal equations is quite difficult and unfortunately it is consid- 
ered as insignificant. Indeed, if a direct problem is smoothing or damping, 
then after a sufficiently long time period it does not matter if the exact or 
the perturbed model is used. However, the situation is quite different for the 
respective inverse and ill-posed problem for which data and modeling errors 
have a strong impact on the solution. Hence, although our causal diffusion 
model yields similar numerical values as the standard diffusion model (com- 
pare Fig. [3] and Fig. S]), it seems evident that causal diffusion is of practical 
interest for inverse problems related to diffusion. 

Apart from this fact, the modeling and investigation of causal mathemat- 
ical equations is interesting from the pure mathematical point of view. 

The goal of this paper is to investigate to what extent a causal diffusion 
model influences the respective backwards diffusion problem. We show that 
this inverse problem is ill-posed, but not exponentially ill-posed. For this 
purpose we derive basic properties of the causal diffusion model (cf. Sec- 
tion |2] and the appendix) which can be summaries as follows. If v denotes 
the distribution of a substance diffusing with constant speed c and initial 
concentration u, then we have the following analytic representatioiu 

(1) V{;m + S) = (2 7r)-^/2^^(| . |)m^^(| . I ^)^ 

for m G No and s G (0, 1], where v denotes the Fourier transform of v with 
respect to x G M^ and Tn is the solution of 

nW + ^^^^ T^W + T^W =0 t > , 

with initial conditions T7v(0+) = 1 and T^(0+). For example, for A^ = 
1, 2, 3 we have 

(2) Ti(t) = cos(t), T2(t) = Jo(t) and Tsit) = sinc(t) , 

where Jq denotes the Bessel function of first kind and order zero (cf. ap- 
pendix). We note that from property ([1]) and supp(T(-,m -|- s)) = Bm+siO) 



^We will see that causal diffusion is determined by the speed of diffusion c, a time 
period r and the space dimension N. Here we assume c = 1 and t — 1. 



causality can be infered (cf. ([7])). Here T(-,s) denotes the inverse Fourier 
transform of T(| ■ | s). Several important properties of the causal diffusion 
model are derived from these properties which are in strong contrast to the 
standard diffusion model (cf. Section [2] and [3j). 

The respective backwards diffusion problem corresponds to the solution 
of the Fredholm integral equation of the first kind 

Ft{u) = w for given data w, 

where the forward operator is defined by Ft{u) := v{-,T) with f as in ([T]) and 
T > denotes the data acquisition time. We show (for appropriate spaces) 
that the forward operator is injective and that it is compact (cf. Section H]) 

1) if A^ = 2 and T > 2 r and 

2) if A^ > 3 and T > r. 

Here N denotes the space dimension and t the time. We note that the enve- 
lope of the Fourier transform of Ft{u) does not decrease exponentially fast. 
In this sense the inverse problem is not exponentially ill-posed. Furthermore, 
numerical simulations of the backwards diffusion problem are performed (cf. 
Section [5]), which confirm our theoretical results. 

The paper is organized as follows: In Section [2] we present our causal 
model of diffusion and derive those properties of diffusion that are needed 
for this paper. For the convenience of the reader we put the technical part 
that is relevant for Sections [2] and [3] in the appendix. Comparisons between 
the standard diffusion model and our causal diffusion model are performed 
in Section [31 The theoretical and numerical aspects of the backwards diffu- 
sion problem are investigated in Sections H] and O Numerical simulations of 
the inverse problem via the Landweber method are presented at the end of 
Section [51 

2 Causal diffusion and its properties 

We now define causal diffusion for the case of a constant speed c G (0, oo). 
For the more general case we refer to [23] . 

Definition 1. Let c,t & (0, oo), dcr(x') denote the Lebesgue surface measure 
on M^ and [^^^(O)! denote the surface area of the sphere S'r(O). Diffusion 
with a constant speed c is defined by 



(2) / lOij(t) 



v^^^{x,t) = / — — -— -— dcr(x) with Vc,r\t=0 = U, 

J \^R(t)[^)\ 



where Tn{t) '■= n(t)T and 

n{t) G No such that t G (n(t) r, (n(t) + 1) r] , 

and R{t) := c{t — n{t)T). Ifu{x.) = (5(x), then we callGc,T '■= v^t the Green 
function of diffusion. Here (5(x) denotes the delta distributioiu 



on K^ . 



Let c > 0, r > and fc,T(x, t) be defined as in Definition [H It follows 
from induction (cf. Lemma 2 in [25J) that the forward operator 

(4) FT:L\R^)-^L\R^),u^v^^r{^,T) (T > fixed) 

is well-defined and that ||m||li = HfcrlUi) i-e. causal diffusion satisfies the 
conservation law of mass. In order to analyse the properties of the forward 
operator in Section HJ we derive the Fourier representation of the Green 
function Gc^r of causal diffusion. In this paper /(k) and J-'{/}(k) denote the 
Fourier transform of x G M^ H- /(x). Our definition of the Fourier transform 
and the respective Convolution Theorem are formulated at the beginning of 
the Appendix. 

Remark 1. The reader may object that the above definition of causal diffu- 
sion is not derived from first principles and that it does not look very sim- 
ilar to standard diffusion. Because of property ||m||li = II'^c.tUli, it follows 
that Vc^T satisfies a continuity equation and, as shown in the introduction 
of l25f . Vc^T satisfies approximately Pick's law. Hence it is reasonable that 
the causal diffusion model follows from first principles under the side condi- 
tion of causality. The derivation of our model from microscopic equations is 
intended to be carried out in the future. 

Moreover, because standard diffusion satisfies a strongly continuous semi- 
group property with respect to time and, as shown in Theorem [H below, our 
causal diffusion model satisfies a discrete semigroup property with respect 
to time, it is evident that both models yield similar numerical results under 
appropriate conditions (compare Fig. and Fig. ^. 



Theorem 1. Let G^^t o,nd T^ be defined as in DefinitionlJ\ and (f^l) (cf. 
Appendix) , respectively. Then 



TAr(|k|cS 



(5) G,,.(k, ,) = -iliL_J_Z for keW\seiO 



T 



and fis{A) := /^Gc,r(x, s) dx defines a positive measure on the Borel sets. 
Moreover, ^ and [2D\) (cf. Appendix) hold. 



^Our notation of the delta distribution is specified at the beginning of the Appendix. 



Proof. We note that s G (0, r] implies n{s) = and R{s) = ct. Moreover, 
we have 

I /(x')da(x')= I /(x + i?y)i?^-Ma(y), 

5fl(x) 5i(0) 

|5R(0)| = |5i(0)|i?^-^ 
From these facts and Definition [1] with m(x) = S{x.), it follows that 

'''^''"•''= I ^^w'^'^'^y K(o)i ^^<y'' 

5fl(s)(x) 5i(0) 

//s(^) is a positive measure, since Gc,t{'^, s) is a positive distribution. 

To determine the Fourier transform of Gcr, we use the following series 
representation derived in Lemma [2] (cf. Appendix): 






(|k|cs 



>2i 



for s G (0, r] with ao = 1 and 



""^^ ~ (2j)! AT . (iV + 2) • (iV + 4) • • ■ (iV + 2 J - 2) ^o^^^I^- 

In Theorem it is shown that Ti(t) = cos(t), T2(t) = Jo(|k|cs) and fl2^ 
holds. T3(t) = sinc(t) can be concluded from Theorem [HI too, or alternatively 

from 

1 1 ■3-5---(2j-l) _ 1 

""^^ ~ J2])\ 3-5-7---(2j + l) " (2j + l)! " 
This concludes the proof. D 

The following theorem together with Theorem [T] provides us with a com- 
plete description of causal diffusion and a mean to compare causal and stan- 
dard diffusion in the k — t— space (cf. Subsection 13. ip . 

Theorem 2. Let Vc^r, u o-nd Gc^r be defined as in DefinitionUl Moreover, 
let Sr denote the space convolution operator with kernel G{'x, t), i.e. SrU := 
G{-,t) *x u for every u G L^(R"). Then we have 

Vc,r{-,t) = S!^ SsU for t = Tm + S,Se{0,T], 

which is equivalent to 

(6) V(-,t) = (27r)-^/2T^(| ■ Icrrr^d ■\cs)u. 



Proof. The claim follows by induction. Let m = 0. Then t = s G (0, r] and 
from Theorem [H we get 

(a,.(-,.)*xu)(x) = y' J ^^^l^^^^^x-x')dx'da(y) 

RN Si{0) 
Si(0) 

Now we assume the induction assumption 

v{-,r + Trn^^) = S!;-'SrU ioT r G (0, t] . 

Let t = s + Tm with s G (0, r]. From the induction assumption with r := r, 
Theorem [1] and Definition [H we infer 

S!; S, «(x) = 5, 5™ «(x) = 5, i;(x, r„) 

/■ /■ (5(x' + csy) , \A 'A / \ 

= J J — \SlO)\ — '"(^-^^^rn)d^da{y) 

R^Si{0) 






Si(0) 

since n{t) = m. 

Finally, the representation (jH]) follows from Theorem [T] and the Convolu- 
tion Theorem (12T|) (cf. Appendix). This concludes the proof. D 

Remark 2. According to Theorem\^ the family of operators {S^-^ | ?7i G No} 
is a discrete semigroup, i.e. So is the identity and 

Sr^+r„ = Sr,„ Sr„ for m, 71 G Nq . 

Formally, the limit r — )■ leads to a continuous semigroup {St\t > 0}. 
Indeed, in Subsection \3.2\ we show for the case N = 2 that the limit r — )■ 
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under the side condition |-^ = const, yields the standard diffusion model and 
thus limT-_j.o{S'T-„ | m G Nq} is a strongly continuous semigroup on [0, oo). A 
consequence of this limit process is that c -^ oo, i.e. this limit diffusion 
process is not causal. Because of these facts, we denote the Green function 
of standard diffusion by Goo,o- 

Corollary 1. The diffusion model defined by Definition^ satisfies causality 
condition 

(7) supp{Gc,r) C {(x,t) G M^ X [0,oo) I |x| < ct} . 



Proof. We recall that the space convolution of two distributions with compact 
support A and B is well-defined and its support lies in A+B := UxeA{x} + i?. 
Moreover, we note that 

Br,{0) + BrM = Br,+rM, 

where Br{0) denotes the open ball of radius R and center 0. 

Let t = s + Tm with m & Nq and s G (0, r]. From Definition [1], it follows 
that supp(G'c,t(x, r)) C Bcr{0) for r G (0, r] and thus by Theorem [2] the 
support of Gc^ri^, t) lies in Bcs+cmriO) = -Bct(O), which proves the claim. D 

The following corollary provides us with an alternative definition of causal 
diffusion and it can be used to compare causal and standard diffusion in the 
X — t— space (cf. Subsection I3.2p . 

Corollary 2. Let Tm for m G No and Vc,t be defined as in DefinitionUl with 
n G L^(M^). The function 



^c,T\ 1 'm 



Tm + s) for s G (0, r] 
solves the wave equation 

(8) _ + ^_^__c2V='t; = on (0,r] 



with initial condition^ 

dv 
(9) ^^(■,0)=t;e,.(-,rj and — (.,0+) = 0. 

Here fc,r('5 0) (m = 0) is understood as the initial distribution u E L 
Proof. From Corollary [5] (cf. Appendix), it follows that 

a(|k|cs)^ |k|cs c/(|k|cs) 

where k and c are fixed. From this, and 

i)(k,s) = T^(|k|cr)"^T^(|k|cs)u(k), 
we obtain 

OH {N-l)dv 2-2,1,2^ n 

as^ s OS 



ltTU>N\ 



''it can be shown that g"'^ (■, Tm— ) ^ = — fe^(-,Tm+), i.e. Vc.t is not continuous 
at the set of time instants {Tm \ in G N} where the semigroup property holds. This is in 
strong contrast to standard diffusion. 



But this is equivalent to equation (|H]). The first initial condition follows from 
Tjv(0+) = 1 (cf. Corollary E]) and Vc,riKr^m+) = ^(k) T^(|k| ct)"^: 

v{k, 0+) = u{k) T^(|k| cr)'" T^(|k| cO+) = {;,,.(k, r^+) . 

Finally, the second initial condition follows from T'(0+) = (cf. Corollary E]): 

|^(k,0+) = ^(k)T^(|k| err ^(|k|c.)|,=o+ = 0. 
This concludes the proof. D 

3 Standard diffusion versus causal diffusion 

In the following we compare standard and causal diffusion in the k — t— space 
and the x—t— space. As explained in Remark[21 we denote the Green function 
of standard diffusion by Goo.o, since c = oo and r = 0. Similarly we use the 
notation 

foo,o := Goofi *x u for u e L^(M^) . 

3.1 Comparison in the k — t— space 

First we recall the definition of the Green function Goo,o of standard diffusion 
and establish the link relation 



(10) Dc 



2N 



between the diffusivity Dq of standard diffusion and the parameters c and r 
of causal diffusion. Then we compare the Green function of both processes 
in the k — t— space. 

It is well-known that the Green function of standard diffusion reads as 

follows (cf. e.g. nziiiaiaiiiiiES]) 

G^,o(x,t) := {An Dot)-''/' exp ("J^) (x G M^, t > 0) 

and that its Fourier transform with respect to x is given by 

(11) Goo,o{Kt) = i2n)-^/' exp{-Do\k\H) (k G M^, t > 0) . 

For sufficiently small |k| (and t = r) we have 

G«,,o(k,r)^(2 7r)-^/^ (l - Do |k| V) . 
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For causal diffusion we get a similar approximation from Theorem [T] together 
with Lemma [2] (cf. Appendix), namely 



2N 

Comparison of these first order approximations yields the link relation (ITOj) . 
Remark 3. Because of 

2a^gaiog2 j^^ ae^, 

the function 



(12) G#(x,t):=(4 7rDot)-^/='2V ^^o*j (x G R^, t > 0) . 

satisfies the standard diffusion equation with diffusion constant D^ := Dq/ log(2), 
i.e. 

(13) G#(k, t) = (2 7r)-^/2 g^p ^_jj^ |]^|2 ^^ _ 

VKe use this (perturbed) diffusion model as a third reference model. 

For the rest of this subsection we focus on the case A^ = 3 for which we 
have (cf . Theorems [1] and [2]) : 

(14) G,,r{k,T^ + s) = {2 7Ty^/^smc'^{\k\cT)smc{\k\cs) sG(0,r]. 

We see that the function k \-^ G'c,r(k, t) (t > 0) is not C°°, since the necessary 
condition 

3C>0VmGNVkGM^: |Gc,r(k,t)| <C(l + |k|)-™ 

does not hold (cf. Paley- Wiener-Schwartz Theorem in [2D]). Moreover, it is 
easy to see from flT^ that 

1) k I— )■ G'c,r(k, t) has a discrete and infinite set of zeros, 

2) k I— 7- Gc,T(k, ti) and k i— )■ G'c,T(k, t2) have the same zeros if and only if 

(ti - tsi/r G No, 

3) some of the zeros move with speed c during the time intervals (r^-i, t^] 
{m G N) such that at the time instants t = Tm~i and t = Tm the same 
set Z of zeros occur. The set Z is given by {k G M'^ | G'c,r(k, r) = 0}. 

9 



This behaviour is in contrast to standard diffusion, since k i— )■ Goo.olk, t) is 
C°° and has no zeros at alL However, since Gc^ri'^t) has compact support 
for each t > 0, the Paley- Wiener-Schwartz Theorem imphes a behaviour of 
such type for causal diffusion. 

Now we perform a numerical comparison. 

Example 1. Let N = 3, c=1,t = 1 and Dq he defined as in / fiOj) . For 
these parameters Fig. [I] shows a numerical comparison of the Green function 
of causal diffusion [I4\ ) with the Green functions / flT]) and l[TS\) of standard 
diffusion. This and further numerical experiments indicate that 

i) h h^ G#(k,t) is closer to k i— )■ Goo.olk, t) than k i— t- Gc,T(k, t), 

a) if t is large, then (^^(k, t), Goo.olk, t) and Gcrlk, t) are very small for 
large |k|, 

Hi) k I— 7- Gcrlk, t) has a discrete set of zeros, in particular it is not mono- 
tone. If n is even, then k i— )■ Gcrlk, t) oscillates around zero. 

This behavior indicates that modeling errors are a serious issue for the back- 
wards diffusion problem. 

3.2 Comparison in the x — t— space 

In the following we demonstrate that for an appropiate parameter set a dis- 
cretization of the standard diffusion equation can yield a similar results as 
the causal diffusion model introduced in Definition [1] 

In order to keep the following formulas and equations short, we focus on 
the two dimensional case. Consider the diffusion of an image with size of 
pixel (Ax)^ and size of time step At := r. We use the notion 

v^j := v{i Ax, j Ax, Tm) for i,jEZ and m G Nq. 

If the length of an image pixel Ax satisfies (cf. Definition [T]) 

R{t) = ct = Ax , 

then we can use the (rough) approximation 

/(xO 



|x-y|=i?(r) 



d(riy) ~ {fi+i,j + /i-ij + fij+i + /i,j-i)/4 . 
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With this discretization the causal diffusion model (E]) is equivalent to 



(15) 






«j 



r 



17 



^i+i,j 






+ vl 



-i,i 



Ax2 



+ 



-'i,i+l 



2<, 



+ w 



«j- 



Ax2 



which is the Forward Euler method of the classical diffusion equation. The 
classical diffusion equation can be obtained for r — )■ under the side condi- 
tion Aa;^/(4r) = const, i.e. the diffusivity corresponds to Dq = Ax'^/{2 Nr) 
with N = 2. Carring out this limit process yields c = cxd, i.e. the diffusion 
speed can be interpreted as infinite. In particular, this shows that the dis- 
crete semigroup {Sr^ \m G Nq} (cf. Theorem [2]) converges to the strongly 
continuous semigroup of the standard diffusion equation. 

Remark 4. Similarly, the discretization of the wave equation ^ for N = 2 
yields the Forward Euler method of the classical diffusion equation, since 



n+l 



dfv H — dtv = -^ 
r 



-M,+ 



in-l 



+ 



1,3 



/n-1 



T^ 



dtv 

T 



Here the CFL condition is satisfied for the discretization At := r and Ax : = 
CT (cf. 15]). However, to obtain the Forward Euler method we have to ne- 
glected the second condition in ^, i.e. 



dv 

a' 



,Tm + ) 







for all m G N . 



The following numerical example indicates that for sufficiently large time 
t the forward Euler method (with fine space discretization) can be considered 
as a noncausal approximation of the causal diffusion model. 



^m/s, R := 10 ^m and r := R/c. 
To the parameters c and r (with 



Example 2. Let N = 2, c = 6.3 ■ 10 
Here we use the notation R := -R(t) 
N = 2) of causal diffusion corresponds the diffusion constant Dq := c^ r/4 = 
1.575 ■10~^m^/s of standard diffusion. The initial mass distribution is shown 
in Fig. 0. To calculate Vc,t defined as in Definition Ui the circles with radius 
R were discretized by 65 points (cf. Fig. \^. The noncausal distribution 
f 00,0 'W^'^s calculated via the forward Euler method for the standard diffusion 
equation. To guaranteed the convergence of this scheme, the discretization 
was chosen as Ax := R/10 and At := -oWrF ^'^^h that 



27VAt 



> D. 
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holds. A time sequence of Vc^ri'^t) and Voofi{-,t) for the time instants t = 
|, ^, r, ..., 2r is visualized in Fig. and Fig. [7} respectively. As ex- 
pected each distribution Voo,o{-,t) is very smooth, in contrast to the distri- 
bution Vc,T{-,t) of causal diffusion. No edge or corner appears in the case 
of standard diffusion. Although it is not visible, in contrast to Vc,T{-,t), the 
support ofVoQfi{-,t) does not lie within the image. For this example, Voofi{-,t) 
and Vc^ri'jt) are very similar after a time period of about t = 3r. Again, we 
note that this means that modeling errors for the backward diffusion problem 
is an issue. 

4 Basic properties of the forward operator 

The calculation of a diffusing substance over the time period T with initial 
concentration u G L^(]R^) corresponds to the evaluation of the forward op- 
erator dl]). We define this as the direct problem and consider the estimation 
of the initial concentration u from appropriate data w. That is to say the 
solution of the Fredholm integral equation of the first kind 

(16) Ft{u) = w for given data w. 

This inverse problem requires the knowledge of c, r and T. In this section 
we investigate the properties of the forward operator and in the subsequent 
section we discuss and perform numerical simulations of the inverse problem. 
We use the notation: 

Definition 2. a) Let T > and Qq be an open subset of Br{0) (r > 0). Then 
we define Qt '■= -Sr+cT(0). Here Br{0) denotes the open ball with center 
and radius r. 
b) L^(M^) is defined as the space of L"^— functions with compact support in 

R^. 

Theorem 3. Let T > and Gc,r be as in DefinitionUl The sets of zeros of 
Gc^r{',T) is discrete and countably infinite, and the operator Ft : Ll{M.'^) — )■ 
Lj(]R^) is injective. 

Proof, a) For t = r„ + s with n G Nq and s G (0, r], we have Gc,t(M, ^n + s) = 
Gc,t{^,t)"' Gc,t{-,s). Hence it is sufficient to show that the sets of zeros 
of Gc,t{-,s) is discrete and countably infinite. Assume that the function 
/ : k I— )■ Gc^t{',s) vanishes on a non-empty M with an accumulation point 
ko G M. Because / has compact support, it can be extended to an analytic 
function f^^t : C^ ^ C^ such that f^^ti^) = /(k) for k G M^ (Paley- Wiener 
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Theorem). Since fextO^) = for k G M and ko G M is an accumulation 
point, /exi is the zero function. Thus M must be discrete. (This fact can 
also be concluded from Theorem [91) That the set of zeros of Gc,t{-,s) is 
countably infinite follows from Theorem |9] and the fact that Ti(s) = cos(s) 
and T2(s) = Jo(s) have countably infinite zeros. 

b) For the injectivity of Ft- Since u and Gc^r{-, T) have compact support, u 
and Gc^r{',T) exist and the Convolution Theorem holds (cf. Theorem 7.1.15 
in |20]). Hence 

T{FTiu)} = G,,riKT)u 



which implies that 
is equivalent to 



Ft(u) = ^ u = 



Gc,ri-,T)u = => u = 0. 



From G^ri', T) u = and part a) of the proof we infer that u vanishes on a 
non-empty open set M. Because u has compact support, the Paley- Wiener 
Theorem implies that u can be extended to an analytic function Uext on C^ 
satisfying {iexi(k) = for k G M. Therefore Uext is the zero function and 
consequently u vanishes. This proves that Ft is injective. D 

Theorem 4. The operator Ft : L^(M^) — )■ L^(M^) is positive, linear and 
self-adjoint. 

Proof. First we show that Ft : Ll{R^) -^ Ll{R^) is well-defined. Because 
L^(R^) is a subspace of a Hilbert space, it is a Hilbert space, too. If m G 
L2(R^), then u G Ll{R^) and thus Ft{u) G L\R^). Since Gc,r and u have 
compact support, their convolution exist and it has compact support (cf. 
Theorem 7.1.15 in [20]). Hence we obtain Ft(L2(R^)) C LJ(R^). According 
to Parseval's formula and 

I (2 7r)^/2 G(k, t)\<G for some constant G , 

which follows from Theorem [9] (cf. Appendix) and Theorem [H we have 

\\Ft{u)\\1. = {2nf\\J^{FT{u)}\\l, = {2nrG J\G{Kt)u{k)\'dk 

— C* 11^11^2 < OO , 

i.e. Ft{u) G L^(R^). Hence the operator is well-defined. 

The positivity and linearity of the operator Ft follows at once from Def- 
inition [1] and Theorem [21 respectively. 
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Since Gc,r{^ — x', T) = Gc,r{'^' — x, T), it follows that 
{Ft{u),w)l'2= [ /"a,r(x-x',T)M(x')u'(x)dx'dx=(M,F:rH)L2 



for w G Lj(M. ) and thus Ft is self-adjoint. This concludes the proof. D 

Theorem 5. If N = 1 and T > 0, then Gc,t{'iT) is a discrete and positive 
measure and the operator Ft '■ L?'{VLq) — > L'^{VLt) is not compact. 

Proof. Without loss of generality we set c = 1 . According to Theorem [1] we 
have for s G (0, r]: 

Gc,t{x., s) = J-'~^{cos{k s)}{x) = - [S{x — s) + 5{x + s)] , 

which implies that Gc,r{x,T) is a convolution of positive distributions with 
singular support. Therefore Gc,t{-,T) corresponds to a discrete and positive 
measure. That Ft is not compact follows from the fact that 



5 



FT = {Rr + LrT{,Rs + Ls) for T = T^ + s 
where Rg, Lg : L^(]R) — > L^(R) are noncompact operators defined by 
Rs{u) := u{- — s) and Ls{u) := u{- + s) . 

D 

Theorem 6. If N = 2 and T > 2t, then G^,r{-,T) G ^^(R^) and the 
operator Ft '■ L'^{VLq) — > I?'{VLt) is compact. 

Proof. Without loss of generality we set c = 1. Let T = Tm + s with m > 2 
and s G (0, r]. From Theorem [T] together with | JqI'")! < 1 and the asymptotic 
behaviour (12^ of Jq (cf. appendix), we get for N = 2: 

oo 

\\G,A-,T)\\hiR^) = '2n I Jl'-{rT)Jl{rs)rdr 



oo 

< A-\ / — dr, 



where 



ITT s J r^ 

M 



M 

A:=2tt / j2™(rr)J^(rs)rdr <cx) 
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and M > 0. Because of J^ 1/r^ dr = 1/M, we arrive at 

2^ 1 

i.e. Gc,r{-,T) lies in L2(M2). Consequently, Gc,r(-,T) lies in ^^(R^). The 
compactness of the operator Fj- (for N = 2 and T > 2 r) follows from 
Theorem 8.15 in [5]. D 

Theorem 7. //A^ >3andT>T, then Gc,r{-,T) E Ll{R^) and the operator 
Ft : L'^i^o) — )■ L'^iS^r) is compact. 

Proof. Without loss of generality we set c = 1 . Let T = Tm + s with m G N 
and s G (0,r]. From Theorem [1] and the estimation fl27|) in Theorem [9l it 
follows for A^ > 3: 

oo 

lia,.(-,T)||i.(^.) = \Sm\ lrJr{rs)T%{rs)r''-'dr 



oo 



^5)(m+l)(7V-l) 
M 

where 

A:= |5i(0)| /"T^'"(rs)T^(rs)r^-Mr<oo 



and M > is sufficiently large. Because of J^r~*^^~^Mr = —jj;^^M~^^~'^\ 
we end up with 

/ ^2(m+l) \ 

i.e. Gc,r(-,T) lies in L'^{R'^). As a consequence, Gc,r(-,T) lies in L2(M2). 
The compactness of the operator Ft (for N = 2 and T > 2 r) follows from 
Theorem 8.15 in [3]. D 

From the Paley- Wiener-Schwartz Theorem (cf. [20]), it follows that the 
Moore-Penrose inverse F^ is uniquely defined by 

F^:={FlrFJ ^T = T^ + s,se{0,T]) 
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with 

(17) ^{Fj(«;)}(k) := ^^xno(s)(k) for k G Qo{s) . 

G(k, s) 



Therefore if the data hes in 

n{FT) ■.= {we l2(m^) 



"^ e L^rM^^ 



G(-,r)-G(-,.) 



then the initial concentration u can be estimated in principle. In contrast 
to standard diffusion G(k, s) has countably infinite and discrete zeros (cf. 
Theorem [3]) • Hence it follows: 

Corollary 3. A necessary condition for w G TZ{Ft) is that w has a zero of 
order > m at k^, if Gc^t{',T) has a zero of order m at k^,. 

According to Theorem IH] for A^ G N we have 

T^(t) X t(^-i)/2 f^^ ^ ^ ^ 

and thus the envelope of k i-t- G^rlk, T) decreases as 

k^ar|k|(L^/-J+i)(^-i)/2 for |k|^oo, 

where 

a^ := (cr)L^/-J(^'i)/^ (c (T - [^^/rj r)^''-'^' . 

Here \_a\ denotes the largest integer < a (and m = [T/rJ , s = T — t [T/rJ). 
Hence we get: 

Corollary 4. If N = 2 and T>2TorN>3 and T > t, then the inverse 
problem i[T5\) is ill-posed, but not exponentially ill-posed. 

We end this section with a remark about the technique of time reversal. 

Remark 5. For the special case T G (0, t], it follows from Corollary {^ that 
the inverse problem l[T5\) depends continuously on the data if the additional 
data W2 '.= — |^(-,T) is known. More precisely, the solution can be calcu- 
lated by 

^ /, N /, N w(k) ,, , Woi^) 



T^(|k| cT) -""^ ^|k|cT'^(|k|cT)' 

where xa is defined as in / f^) and {A, B} is a covering ofW^ such that AO 
i? = and Tat (I ■ | T) and T^(| -IT) do not vanish on A and B, respectively. 
Here we have used the fact that the zeros of T ^ and T^ are of order one 
which follows from CorollarylE and Theorem\^in the appendix. 
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5 Simulation of the inverse problem 

5.1 Simulation of data via a particle method 

In order to avoid an inverse crime we calculate the synthetic data for the 
inverse problem by a particle method (cf. [16j). One of the advantages of 
a particle method (as long as no mass flows over the boundary) is that the 
total mass is conserved. For simplicity we focus on the 2D— case and drop 
the subscripts c and r in G^r and Vc^r- 

The particle method 

The initial distribution u is approximated by an image, i.e. a piecewise 
constant function with quadratic pixels of length Ax. At time instant t = 
Tn-i {n G N) the mass concentrated in a pixel separates in M parts and each 
part propagates on a stright line with constant speed c in a randomly chosen 
direction d during the time period r. Here the directions are chosen with 
equal probability out of the set 

{A{^) ei I y; = 0, vr/M, . . . , (M - 1) tt/M} , 

where ei := (1,0)^ and A[ip) denotes the matrix that rotates the argument 
about the angle if in positive direction. This kind of data simulation allows 
that more than one "particle" go in the same direction such that a special 
type of noise is included in the simulated data. To each image pixel is then 
associated the number of all particles that lie within the pixel multiplied by 
1/M. 

Noise 

In order to avoid an inverse crime we perturbed the length of the radius 
R{t) = cr by ±0.25% of its original length (uniformly distributed perturba- 
tion). In addition, uniformly distributed L^— noise with positive mean value 
were added to the simulated data. As noise level we have chosen 5 = 0.005 
(0.5%). 

Convergence of the particle method 

In the following we denote by G'[M](x, t) the simulated distribution with 
initial distribution 



5[M](x) : = 



1 if max(|x|, \y\) < ^ 
elsewhere 
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From analysis it is known that 

(18) 5[M](x)''^5(x) and G[M]{^,t) ^'^ G{^,t) in V\M.^) . 

Here G denotes the Green function of causal diffusion (cf. Definition [1]) and 
P'(M^) denotes the space of distributions on M^. We now show that the 
algorithm described above for an initial distribution u provides us with an 
approximate solution of Ft{u), were F^ denotes the forward operator (jlj). 

Theorem 8. Let u E Ll(M?) and v{-, t) := G{-, t) *^ u. For 

v[M]{;t):=G[M]{;t)*^u (t>0), 

it follows that 

v[M]{-,t)^-^v{;t) m L\R^). 

Proof. Since the space G^{M.'^) is dense in L^(M^), we assume without loss 
of generality that u G C^(M^). We have 

||^[M](-,t)-i;(-,t)|Ui = y"|/[M](x)|dx 

with 

(19) /[M](x) = y"[G[M](x',t)-G(x',t)]7i(x-x')dx'. 

R2 

The function f[M] is an element of ^^^^(M^), since G[M](-,t) - G{-,t) has 
compact support and u G C^(R^) (cf. Proposition 32.1.1 in [I5]). From 
(IT8|) together with u G C^(M^), it follows that the right hand side of (p|) 
converges pointwise and uniformly to zero on compact sets. This together 
with the fact that f[M] has compact support implies 

\\v[M]{.,t)-v{;t)U^'^0. 

As was to be shown. D 

5.2 Numerical solution of the backwards diffusion prob- 
lem 

For solving the inverse problem we use the Landweber method (cf. e.g. [TU | [2 ^ 
[221 ESI). Since Ft : Ll{R^) -^ ^^(M^) is a positive, linear and self-adjoint 
operator (cf. Theorem Hj) the Landweber method reads as followqj 



^For simplicity, we write u„, w^ instead of u„[Af], w^[M]. 
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where u denotes the relaxation parameter, w^ denotes the noisy data and P 
denotes the orthogonal projection onto 

7^(P) = {u G L^ I n > 0} . 

The use of the projection operator guarantees that the solution is a posi- 
tive (mass) distribution. As parameter choice rule we use the discrepancy 
principle, i.e. the iteration is stopped as soon as 

\\FT{Un+l)-w'\\L2<Tl6 {r]>2) 

is true. The relaxation parameter was chosen as 

|2 

Il2 



Ul 



1 WFtM-w^ ^ 
4\\Ft[FtM-w^]\\1, 



In order to avoid an inverse crime, the data w^ is calculated by the par- 
ticle method (M = 65) described above and the calculation of the Forward 
operator F^ in each iteration step is performed by integrals over circles. Each 
circle is discretized by 50 points. 

We now present two simulations of the backwards diffusion problem for 
T = T and T = 3 r, respectively. 

Example 3. Consider the initial distribution shown in Fig. O This image 
consists o/ 682^ quadratic pixels of length Ax := 1/681. As characteristic 
parameters of causal diffusion we have chosen c = 1 and r = 8 Ax/c. Hence 
the characteristic radius R{t) is 8 times Ax. As described above, the length 
of the radius R{t) was randomly perturbed by ±0.25% of its original length. 
The data acquisition is performed at time T = t and 0.5% (uniformly dis- 
tributed) L"^ — noise was added to the simulated data. The numerical results 
are visualized in Fig. and Fig. 0. ^45 expected, the estimation of large 
structures is much better than for smaller ones. Since the data acquisition is 
performed at a quite early time the estimation works wellu The Discrepancy 
principle stops optimally for rj = 9.4 after 6 steps. 

Example 4. We consider the inverse problem from Example again, but 
for the later data acquisition time T = St. The Discrepancy principle stops 
optimally for rj = 5.9 after 6 steps. For this situation the forward operator 
is compact. As Fig. [^ shows it is not possible to restore the edges of the 
question mark, since the data are to much "smooth". This result reflects the 
ill-posedness of the problem. 



^Cf. Remark [i 
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6 Appendix 

The delta distribution 

We use the following notation for the delta distributions. Let N & N and 
X G M^. Then 6{x) satisfies 

/(x)5(x-Xo)dx = /(xo) for /gC,(M^). 

r 

Here Cc(M^) denotes the set of continuous funtions with compact support. 
Since S{x) has compact support, Cc(M^) can be replaced by C(M^). In this 
notation the dirac measure /i^ (cf. [2Z]) reads as follows 

where 

(20) Xa(x) denotes the characteristic function of the set A C 

In case A^ = 1 we use the notation 6{x) instead of 5(x). 

The Fourier transform 

We use the following notation for the Fourier transformation: 

/(k) :=-F{/}(k) := (2 7r)-^/2 ^ e^'^V(x) dx 
^(x) :=^-i{^}(x) := (2 7r)-^/2 [ e-'^"-- g (k) dk 



oN 



for f,gE L}{EJ^). Here k G M^ is called the wave vector. In this notation 
the convolution theorem reads as follows 

(21) :F{h]:F{h] = (2 7r)-^/2^{A *. f,} A, h G Li(M^) . 

Special functions 

We define the function sine : M — )■ M as the continuous extension of x G 
R\{0} I—)- sin{x)/x and recall that the Bessel function of first kind and order 
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zero has the series representation (cf. [12]) 

(22) jo(x) = J2i-^y 



X- 



J 



2 



i=o ^ -^ 

In order to derive an analytic representation of the Fourier transform of 
the Green function of causal diffusion, we need the following two lemmata. 

Lemma 1. For N eN with N > 1 and j E N, let 

7r/2 

IiN,j) := I sin-'((/?) cos^'^^^) ^^ _ 

-7r/2 

If i is odd, then I{N,j) = and if j is even, then 
HN,j) l-3-5---(j-l) 



(23) 



IiN,0) N-{N + 2)-{N + A)---{N + j-2) 



Proof. If J is odd, then sm^{(p) cos'^~^{ip) is an odd function and thus I{N,j) 
vanishes. 

Now let j be even. We perform a proof by induction. 

i) Let j = 2. Integration by Parts yields 

/(iV,2) = l Tcos"-(^)d»,-'<''-°' 



N J ^'^' "^ N 

-it/2 

ii) We now assume the induction assumption for j = m — 2, i.e. 

1 . ^ . r^ . . . (m — "{] 
^W'"^^'- iV.(iV + 2).(iV + 4)...(A,n-4) ^'^-°' 

and prove (123]1 for j = m. Integration by Parts yields 

771 — 1 

I(N, m) = UN, m-2). 

Employing the induction assumption to this result leads to (l23l) with 
j = m. This concludes the proof. 

D 
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Lemma 2. Let 

(24) 



'^N{t):=J2(-^y-('^rt^' for tG(0, 

j=0 



oo 



with Co = 1 and 
1 

02. 



l-3-5---(2j-i: 



'^ (2 j)! A^ ■ (A^ + 2) ■ (iV + 4) ■ ■ ■ (AT + 2 j - 2) 
r/ie series [24^ is absolutely convergent and 



U e N) . 



5(x + csy) 



|5i(0) 



da(y) 



(2 7r)^/2 



for xGM^\sG(0 



r 



5i(0) 



Proof. That the series representation fl24)) converges absolutely follows at 
once from the Quotient Criterion. 
Let X, k G M^. From 



J-{5(x + aiy)}(k) = (2 7r)-^/^e 



-Af/2„-iai(k-y) 



[tti > constant) 



and 



f e^'^^y^daiy) = f e"^ l'^l("^-y)da(y) (as G C constant), 



5i(0) 

it follows that 
g(k,s) := J" 



Si(0) 



S{- + csy) 

\sim 



My) } (k) 



[SiW 



Si(0) 



g-i(ei-y)|k|cs 

(2 7r)^/2|5i(0)| 



da(y) 



for s G [0, r] and k G M^. Instead of ei we can also use anyone in {e2, 63 , . . . , ejv}- 
Expanding the exponential function yields 



(25) 
with 



(2vr)^/^^(k,.) = 5^(-l)^i^L^d. 



j=0 



(2j)! 



2j 



do 



(ei ■ y) 



d(7(y) for j eN. 



' ■ ' |5i(0)| 

5i(0) 

We see at once that dj = if j is odd and do = I. For the convenience of the 
reader, we consider the cases A^ = 1 and A^ > 1 separately. 
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a) For A^ = 1 we have ei = 1 

J da{y)^ J{5{y-l) + 5iy + l))dy and |5i(0)|=2, 



5i(0) 

and thus 



,_!-' + (—1)-' _ J 1 if J is even 
^ ~ 2 ~ 1 if j is odd 



Inserting this into the series representation yields 

A\k\csr^ 



{2nf/mk,s) = J2i-^ 



CO 
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b) Let A^ > 1. For the derivation of the series representation we use the 
following A^— dimensional orthogonal coordinate system (cf. [52] ) 

(r,(^i, ... ,<^7v_i) G [0,oo) X (-7r,7r) x (-7r/2, 7r/2)^-2 

defined by 

xn = r sin(v?iv-i) , 

Xn-1 = r cos(v?Ar-i) sin(v9Ar_2) , 

Xn-2 = r COs{ipN-l) COs{(^N-2) Sm{ipN-3) , 

Xn-3 = r COs{(pN-l) COs{ipN-2) COS^ipN-s) sin(v5Ar-4) , 

X2 = r cos{ipN-i) ■ ■ ■ cos(v22) sin(v9i) , 
xi = r cos{ipN-i) ■ ■ ■ cos(v92) COs((y9i) , 

with surface measure 

Af-l 



da = r^ ^difi JJcos' \ipi)dipi. 



1=2 



Since e^r ■ e^ = sm{ip]y_i) and 



» N-l 

\Si{0)\= I difi llcos'-\^i)dip^ 



Si(0) '-2 
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we obtain 

N-2 



d,= |d^. n<=-"'fe)d^.^^^||jf 



cos (v9Ar-l)dv9Ar_i 

1=2 



Si(0) 

I{N,0) 

with I{N, j) defined as in Lemmadl From this and Lemmadl we obtain 
d2j-i = 0, do = 1 and 



1 ■3-5---(2j-l) 

N-{N + 2)-{N + 4)---{N + 2j -2) 



d2j = — — fj^j , ^v^ — ^^^ ^ ^^ ^^^ ^ ^ . — ^^ for j > . 



Inserting this into the series (125|) yields the claimed series representa- 
tion. 

This concludes the proof. D 

The following corollary follows form Lemma [2j 



Corollary 5. For N E N. The function T^v defined as in ([^]j satisfies the 
problem 

mt) + ^^^^ T^(t) + T^(t) = t > , 

wt/i initial conditions 

T^(0+) = 1 and T'jv(0+)=0. 
i^ere t = zs a regular singular poinlo of the ordinary differential equation. 

The following theorem enables us to specify the space Fourier transfrom 
of the Green function of causal diffusion for every dimension N and to prove 
some compactness results for the forward operator of causal diffusion. 

Theorem 9. Let N E'N with N > 3 and t > 0. The function Tjv defined as 
in [24\ ) satisfies 

(26) T^(t) = -i^^T^_,(t) 

with 

Ti(t) = cos(t) and T2(t) = Jo{t) . 

6cf. e.g. ng. 
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Here Jq denotes the Bessel function of first kind and order zero. Moreover, 
we have 

(27) |TAr(t)| < C^t-(^-^)/2 for sufficiently large t 

and some constant Cn > 0. 

Proof. The relation between Tat and T'^_2 follows at once from the series 
representation fl2^ . Moreover, 

a) if A^ = 1, then 

1 1 ■3-5---(2j-l) 1 



02, 



(2j)! l-3-5---(2j-l) (2j)! 
and thus T(t) = cos(t) and 
b) if iV = 2, then 

1 1 ■3-5---(2j-l) 1 1 



«2, 



(2j)! 2-4.6---(2j) [2-4-6---(2jP [2i (j!)]2 

which implies T(t) = Jo(t) (cf. f l22|) in the Appendix and [19j). 
In order to proof the estimation we use 

(28) (iV-2)V3 ='^^-2^' 



which follows from the series representation fl24l) . We perform a proof by 
induction. Since cos{t) is bounded and the Bessel function Jo(t) satisfies the 
asymptotic behavior (cf. [1]) 



JoW^y^cos(t-^) for t 



(29) Jo(t) X W — cos(t--j for t ^ oo , 

the estimation holds for A^ = 1 and N = 2. We assume that the estima- 
tion ( 127|) holds and prove 

|Tjv+2(i)| < Civ+2t"^^+^^/^ for sufficiently large t. 
From d^H]) and ([2HD we get 



|T;v+2(t)| 



A^ 



< ^% ^^ |T^(t)| + |T;v-2(t)|) 



Ar(Ar-2)(Cjv + CAr,2) 
^(JV+l)/2 

which proves the claim. D 
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Figure 1: Comparison of |k| i— )• (2 7r)^/^ (^^(k, t) (dotted red line), |k| i— ;■ 
(2 7r)3/2G'^o(k,t) (solid green line) and |k| ^ (2 7r)3/2 G^^^(k, t) (dashed 
black line) . Here c = 1, r = 1 and t G {r, 9 r}. 
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Figure 2: The left picture shows the initial distribution u and the right 
picture visualizes the discretization of the circle of radius R = R{t) that is 
used to evaluate the circle integrals in Definition [H 
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Figure 3: Visualization of causal diffusion for the time sequence 
(|, ^, r, . . . , 2r) with discretization Ax = 9.6 ■ 10^^ and At = t = R/c. 
The initial distribution is shown in Fig. |2i 
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Figure 4: Visualization of standard diffusion for the time sequence 
(|, ^, r, . . . , 2r) with discretization Ax := -R/10 and At : 
initial distribution is shown in Fig. HJ 
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Figure 5: Numerical solution of the inverse problem with 0.5% uniformly 
distributed L^-noise for data acquisition time T = r and T = 3 r. For T = t 
and T = 3 r, the Discrepancy principle stops optimally for r] = 9.4 (6 steps) 
and 1] = 5.9 (6 steps), respectively. 
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Figure 6: Details to the second row in Fig. [5] (T = r). The left column shows 
the left top part of the big question mark and the right column shows the 
small question mark and the small apostrophe. 
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Figure 7: Details to the third row in Fig. |5] (T = 3 r). The left column shows 
the left top part of the big question mark and the right column shows the 
small question mark and the small apostrophe. 
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